function [imgVolume] = reorder_raw_siemens_mridata(fid,mdh,seq)
%REORDER_RAW_SIEMENS_DATA reorders meas.dat files opened with rdMeas. 
%   Input : reorder_raw_siemens_mridata(fid,mdh)
%   Output: imgVolume(fe,pe,sl,ch)
%
%   fe = frequency encoding
%   pe = phase encoding
%   sl = slice
%   ch = channel
%
%   ***************ATTENTION****************
%   *****Missing lines are zero-filled.*****
%   ****************************************
%
%   ***************ATTENTION****************
%   **Can only process data with one slab.**
%   ****************************************
%
%   ******************ATTENTION********************
%   **Only currently works with simplest GRE data**
%   ***********************************************
%
%   Author: Mahtieu Boudreau
%   Date: March 27th 2013
%   Email: mathieu.boudreau2@mail.mcgill.ca


% Get volume dimensions

numFE = max(mdh.ushSamplesInScan/2);  % Division by 2 because 
                                      % Siemens oversamples their feq. enc.
numPE = max(mdh.ushLine)+1;

%Handle 3D/MultiSlice/MultiPartition data
if (max(mdh.ushPartition) == 0)&&(max(mdh.ushSlice) > 0)
    
    numSL = max(mdh.ushSlice)+1;
    
elseif(max(mdh.ushSlice) == 0)&&(max(mdh.ushPartition) > 0)
    
    numSL = max(mdh.ushPartition)+1;  % Yes I should have renamed this to 
                                      % numPT or something, and put flags 
                                      % later on, but quicker this way.
    
elseif (max(mdh.ushSlice) == 0)&&(max(mdh.ushPartition) == 0)
    
    numSL = 1;
    
else
    
    disp('***Warning, cannot handle max(mdh.ushSlice) && (max(mdh.ushPartition)  both greater than 0***')
    return
    
end

numCH = max(mdh.ulChannelId)+1;

% Pre-allocate imgVolume size.
switch seq
    case 'gre'
        imgVolume = zeros(numFE, numPE, numSL, numCH, 'single');
    case 'mt'
        numMT = max(mdh.ushRepetition)+1;        
        imgVolume = zeros(numFE, numPE, numSL, 1, numCH, 'single');

end

imgVolume = complex(imgVolume,0);

%
numAcqLines = size(fid,2);

switch seq
    case 'gre'
        for ii = 1:numAcqLines
            if ((max(mdh.ushPartition) == 0)&&(max(mdh.ushSlice) > 0))

                imgVolume(:, mdh.ushLine(ii)+1, mdh.ushSlice(ii)+1, mdh.ulChannelId(ii)+1)=fid(:,ii);    

            elseif((max(mdh.ushSlice) == 0)&&(max(mdh.ushPartition) > 0))
                imgVolume(:, mdh.ushLine(ii)+1, mdh.ushPartition(ii)+1, mdh.ulChannelId(ii)+1)=fid(:,ii);  
            
            elseif((max(mdh.ushSlice) == 0)&&(max(mdh.ushPartition) == 0))
                imgVolume(:, mdh.ushLine(ii)+1, 1, mdh.ulChannelId(ii)+1)=fid(:,ii);  
                
            end
        end
    case 'mt'
        for ii = 1:numAcqLines
            ii/numAcqLines*100
            if ((max(mdh.ushPartition) == 0)&&(max(mdh.ushSlice) > 0))

                imgVolume(:, mdh.ushLine(ii)+1, mdh.ushSlice(ii)+1, mdh.ushRepetition(ii)+1, mdh.ulChannelId(ii)+1)=fid(:,ii);    

            elseif((max(mdh.ushSlice) == 0)&&(max(mdh.ushPartition) > 0))
                imgVolume(:, mdh.ushLine(ii)+1, mdh.ushPartition(ii)+1, mdh.ushRepetition(ii)+1, mdh.ulChannelId(ii)+1)=fid(:,ii);  
            
            elseif((max(mdh.ushSlice) == 0)&&(max(mdh.ushPartition) == 0))
                imgVolume(:, mdh.ushLine(ii)+1, 1, mdh.ushRepetition(ii)+1, mdh.ulChannelId(ii)+1)=fid(:,ii);  

            end
        end
end

end

